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i Abstract. We address the classical stellar-atmosphere problem and de- 

scribe onr method of numerical solution in detail. The problem consists 
of the solution of the radiation transfer equation under the constraints of 
hydrostatic, radiative and statistical equilibrium (non-LTE). We employ 
the Accelerated Lambda Iteration (ALI) technique, and use statistical 
methods to construct non-LTE metal- line-blanketed model-atmospheres. 
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On ' 1. Introduction 

o 

, Essentially all our knowledge about the structure and evolution of stars, hence 

about galactic evolution in general, rests on the interpretation of their electro- 
magnetic spectrum. Therefore quantitative analysis of stellar spectra is one of 
the most important tools of modern astrophysics. The formation of the ob- 
served spectrum is confined to the atmosphere, a very thin layer on top of the 
stellar core. Spectral analysis is performed by modeling the temperature- and 
53 . pressure-stratification of the atmosphere and computing synthetic spectra which 

are then compared to observation. Fitting synthetic spectra to the data yields 
basic photospheric parameters: effective temperature, surface gravity, and chem- 
^ I ical composition. Comparison with theoretical evolutionary calculations allows 

. one to derive stellar parameters: mass, radius, and total luminosity. 

The classical stellar-atmosphere problem considers radiation transfer 
through the outermost layers of a star into free space under three assumptions. 
First it is assumed that the atmosphere is in hydrostatic equilibrium, thus, the 
matter which interacts with photons is at rest. Second, the transfer of energy 
through the atmosphere is assumed to be done entirely by photons, i.e. heat- 
conduction and convection are regarded as negligible (radiative equilibrium). 
But the effectiveness of photon-transfer depends on the opacity and emissivity 
of the matter, which are strongly state- and frequency-dependent quantities. 
They depend in detail on the occupation numbers of atomic/ionic levels, which, 
in turn, are determined by the local temperature and electron density, and the 
radiation field, whose nature is non-local. The occupation number of any atomic 
level is determined by a balance among radiative and collisional population and 
de-population processes (statistical equilibrium; our third assumption), i.e. the 
interaction of atoms with other particles and photons. Mathematically, the 
whole problem requires the solution of the radiation transfer equations simulta- 
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neously with the equations for hydrostatic and radiative equihbrium, together 
with the statistical equihbrium (rate equations). 

A stellar atmosphere radiates into circumstellar space and thus is an open 
thermodynamic system; hence it cannot be in thermodynamic equilibrium (TE), 
and cannot be characterized by a single temperature. The idea of "Local Ther- 
modynamic Equilibrium" (LTE) is a hypothesis which assumes that while TE 
does not hold for the atmosphere as a whole, it can be applied in small volume 
elements. In this case the atomic occupation numbers depend only on the local 
electron temperature and electron density via the Saha-Boltzmann equations. 
But this approximation can be valid only in the limit that collision rates dom- 
inate radiative rates, and photon-mean-free-paths are small. Models in which 
the Saha-Boltzmann equations are replaced by the physically more accurate rate 
equations are called non-LTE (or NLTE) models. This designation is somewhat 
imprecise because the velocity distribution of particles is still assumed to be 
Maxwellian, i.e. we can still define a local temperature. NLTE calculations are 
more costly than LTE calculations, however, and it is hard to predict a priori 
whether NLTE effects will be important in a specific problem. Generally, NLTE 
effects are large at high temperatures and low densities, which imply intense 
radiation fields, hence frequent radiative processes (which have a non-local na- 
ture and, in particular, respond to the presence of the open boundary), and 
less-frequent particle collisions, which tend to enforce LTE conditions. 

Abandonment of the LTE assumption leads to a much more difficult model- 
atmosphere problem, because of subtle couplings between the radiation transfer 
and statistical equilibrium equations. The pioneering work by Auer &; Mihalas 
(1969) provided a basic tool for making such models. But the numerical prob- 
lem going from LTE to realistic NLTE models has been solved only recently, 
and is the topic of this paper. We now have much more powerful tools to com- 
pute non-classical models, which account for very complex opacities, and solve 
the radiation transfer equations in more general environments, for example in 
expanding stellar atmospheres. These are the topics of other papers in this 
volume. 

Stellar-atmosphere modeling has made huge progress in recent years. This 
advance was achieved by the development of new numerical techniques for model- 
construction, and by the availability of atomic data for many species. And 
these achievements became possible only with an enormous increase of com- 
puting power. Model atmospheres assuming LTE have been highly refined by 
the inclusion of many more atomic and molecular opacity sources, and effective 
numerical techniques for LTE model-computation have been available for years. 
The progress is most remarkable in the field of NLTE model atmospheres. Re- 
placement of the Saha-Boltzmann equations by atomic rate-equations requires 
radically different numerical-solution techniques, otherwise metal-opacities can- 
not be taken into account at all. Such techniques have been developed with 
great success during the last decade, inspired by important papers by Cannon 
(1973) and Scharmer (1981). The Accelerated Lambda- Iteration Method (ALI) is 
at the heart of this development. Combined with statistical representations of 
line-opacities, we are finally able to compute metal-line-blanketed NLTE models, 
including many millions of spectral lines, with a very high level of sophistication. 



Model Photospheres with Accelerated Lambda Iteration 



3 



In this paper we discuss the basic ideas behind the new numerical methods 
for NLTE modeUng. We begin by presenting the classical model-atmosphere 
problem and introduce the basic equations. Then we focus on the ALI solution- 
method, and our numerical implementation of it. We then briefly describe the 
solution of the NLTE metal-line-blanketing problem. 

2. Statement of the problem and overview of the solution method 

Wc assume plane-parallel geometry, which is well justified for most stars be- 
cause the atmospheres are thin compared to the stellar radius. The only pa- 
rameters needed to characterize uniquely such an atmosphere are the effective 
temperature (Tefj ), a measure for the amount of energy transported through the 
atmosphere per unit area and time (see Eq. 15), the surface gravity (g), and the 
chemical composition. Generalization to spherical symmetry to account for ex- 
tended (static) atmospheres mainly affects the radiation transfer equation and 
is straightforward (Mihalas &: Hummer 1974). A spherical-symmetry variant of 
our code has been written by Nagel etal. (2001). To construct model atmo- 
spheres we mTist solve simultaneously a set of equations that is highly coupled 
and non-linear. Because of the coupling, no equation determines uniquely a 
single quantity - all equations jointly determine each of the state parameters. 
Neverlethess, each of them is usually thought of as determining, more or less, a 
particular quantity. These equations are: 

• The radiation transfer equations which are solved for the angle-integrated 
mean intensities Ji,i = 1, . . . ,NF, on a pre-chosen frequency grid com- 
prising NF points. The formal solution is given by J = AS, where S is the 
source function as defined later (Eq. 11), and A is the transport operator. 
Although A is written as an operator, one may think of it as a process of 
obtaining the mean intensity from the source function. 

• The hydrostatic-equilibrium equation which determines the total particle 
density N. 

• The radiative-equilibrium equation from which the temperature T follows. 

• The particle-conservation equation, determining the electron density Ue- 

• The statistical-equilibrium equations, which are solved for the population 
densities n^, i = 1, . . . , NL of the atomic levels that are allowed to depart 
from LTE (NLTE levels). 

• The equation that defines a fictitious massive-particle density nn, which 
is introduced for convenience in the solution procedure. 

This set of equations has to be solved at each point d of a grid comprising 
ND depth points. Thus we are looking for solution vectors 



ip'a = {ni,...,nNL,ne,T,nH,N,Ji,...,JNF), d=l,...,ND. (1) 
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The Complete Linearization (CL) method (Auer &; Mihalas 1969) solves this 
set by linearizing the equations with respect to all variables. The basic advan- 
tage of the ALI (or "operator splitting") method over CL is that it allows one 
to eliminate, at the outset, the explicit appearance of the mean intensities Jj 
from the solution scheme by expressing these variables by the current, yet-to- 
be determined, occupation numbers and temperature. This is accomplished by 
an iteration procedure which may be written as (suppressing indices indicating 
depth and frequency dependence of variables): 

jn ^ ^*gn + (A - A*)S"-^ (2) 

This means that the actual mean intensity at any iteration step n is computed 
by applying an approximate lambda operator (ALO) A* on the actual source 
function S"" plus a correction term that is computed from quantities known from 
the previous iteration step. This correction term includes the exact lambda 
operator A, which guarantees convergence to the exact solution of the radiation 
transfer problem J = AS. The use of A in Eq. 2 only indicates that a formal 
solution of the transfer equation is performed, but in our application the operator 
is not constructed explicitly. Instead we employ the Feautrier (1964) solution 
scheme (Mihalas 1978) or a short-characteristic method (Olson &; Kunasz 1987) 
to solve the transfer equation as a differential equation. The resulting set of 
equations for the reduced solution vectors 

tp^ = {ni,...,nNL,ne,T,nH,N), d = l,...,ND (3) 

is of course still non-linear. The solution is obtained by linearization and it- 
eration which is performed either with a usual Newton-Raphson iteration or 
by other, much faster methods like the quasi-Newton or Kantorovich variants 
(Dreizler & Werner 1991, Hubeny & Lanz 1992) (see Sect. 5.4.). The first model- 
atmosphere calculations with the ALI method were performed by Werner (1986). 

Another advantage of the ALI method is that explicit depth coupling of 
the solution vectors Eq. 1 through the transfer equation can be avoided if one 
uses diagonal (i.e. local) ALOs. Then the solution vectors Eq. 3 are independent 
from each other and the solution procedure within one iteration step of Eq. 2 is 
much more straightforward. Depth coupling is provided by the correction which 
contains the exact solution of the transfer equation. The hydrostatic equation, 
which also gives explicit depth coupling, may be taken out of the set of equations 
and can - as experience shows - be solved in between two iteration steps of Eq. 2. 
Then full advantage of a local ALO can be taken. The linearized system may 
be written as: 

^, = ^US^d, (4) 

where is the current estimate for the solution vector at depth d and 5if)^ is the 
correction vector to be computed. Using a tri-diagonal A* operator the resulting 
system for is - like in the classical CL scheme - of block tri-diagonal form 
coupling each depth point d to its nearest neighbors d ± 1: 

Id.S'^d-i + f^Jtpd + otdSipd+i = Cd- (5) 

The quantities a, (3, 7 are {NN x A'^A'") matrices where NN is the total number 
of physical variables, i.e., NN = NL + A, and Cd is the residual error in the equa- 
tions. The solution is obtained by the standard Feautrier scheme. With starting 
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values Di = (-cti) and vi = Ci we sweep from the outer boundary of 
the atmosphere inside and calculate at each depth: 

Vd = (/3rf + 7dDd-i)"^(cd-7dVrf_i). (6) 

At the inner boundary we have Djvd = and by sweeping backwards to the 

surface we calculate the correction vectors, first (^i/jjvD = ^ND, and then suc- 
cessively dij^^ = D(^(5'0(^_|_j^ + V(^. As already mentioned, the system Eq. 5 breaks 
into ND independent equations Sil^^ = /3^^C(^ {d = 1, . . . , ND) when a local 
A* operator is used. The additional numerical effort to set up the subdiagonal 
matrices and matrix multiplications in the tri-diagonal case is outweighed by the 
faster global convergence of the ALI cycle, accomplished by the explicit depth 
coupling in the linearization procedure (Werner 1989). 

The principal advantage of the ALI over the CL method becomes clear at 
this point. Each matrix inversion in Eq. 6 requires + operations whereas 
in the CL method (NL + NF + 4)^ operations are needed. Since the number of 
frequency points A^F is much larger than the number of levels NL, the matrix 
inversion in the CL approach is dominated by NF. 

Recent developments deal with the problem that the total number of atomic 
levels tractable in NLTE with the ALI method described so far is restricted to 
the order of 250, from experience with our model atmosphere code PR02. This 
limit is a consequence of the non-linearity of the equations, and in order to 
overcome it, measures must be taken in order to achieve a linear system whose 
numerical solution is much more stable. Such a pre-conditioning procedure was 
first applied in the ALI context by Werner Sz Husfeld (1985) using the core- 
saturation method (Rybicki 1972). More advanced work achieves linearity by 
replacing the A operator with the \I' operator (and by judiciously considering 
some populations as "old" and some as "new" ones within an ALI step) which 
is formally defined by writing: 

Ju = ^uVu , i.e. ^'i. = A^/xi. , (7) 

where the total opacity Xu (as defined in Sect. 3.6.) is calculated from the pre- 
vious ALI cycle. The advantage is that the emissivity r)i, (Sect. 3.6.) is linear 
in the populations, whereas the source function is not. Hence the new op- 
erator ^ gives the solution of the transfer problem by acting on a linear func- 
tion. This idea was proposed by Rybicki & Hummer (1991) who applied it to 
the line-formation problem, i.e. regarding the atmospheric structure as fixed. 
Hauschildt (1993) and Hauschildt & Baron (1999) generalized it to solve the 
full model-atmosphere problem. In addition, splitting the set of statistical equa- 
tions, and solving it separately for each chemical element, means that now many 
hundreds of levels per species are tractable in NLTE (see Dreizler, this volume, 
for our experience). A very robust method and fast variant of the ALI method, 
the ALI/CL hybrid scheme, permits the linearization of the radiation field for 
selected frequencies (Hubeny &; Lanz 1995), but it is not implemented in PR02. 



6 



Werner, Deetjen, Dreizler, Nagel, Ranch, &: Schuh 



3. Basic equations 
3.1. Radiation transfer 

Any numerical method requires a formal solution (i.e. atmospheric structure 
already given) of the radiation transfer problem. We disregard here polarization 
effects. (An implementation of an ALI technique to treat polarized radiation 
is discussed by Deetjen et al. in this volume.) Radiation transfer at any depth 
point is described by the transfer equation: 

±^i^^^ = s,-I,i±^,), Me [0,1], (8) 

written separately for positive and negative (the angle-cosine between the 
direction of propagation and the outward-directed normal to the surface), i.e. 
for the inward and outward directional intensities at frequency u. Here 
is the optical depth (defined in terms of the column mass m used in the other 
structural equations hj dTv = dmxv/ P-, where p is the mass density), and Sj^ is 
the local source function. Introducing the Feautrier intensity 

u,^ = {I,{lj)+I,{-li))/2 (9) 

we obtain the second-order form (Mihalas 1978, p. 151): 

m'^^ = ^^.m-'5.> Me [0,1]. (10) 

We may separate the Thomson emissivity term (scattering from free electrons, 
assumed coherent, with cross-section (jg) from the source function so that 

Sy = Sl + ne(TeJu/Xv^ (H) 

where S'^^ is the ratio of thermal emissivity to total opacity as described in detail 
below (Sect. 3.6.): 5^ = r]y/xv The mean intensity is the angular integral over 
the Feautrier intensity, hence the transfer equation becomes 

IJ- ^ 2 = ~ / ^J'M dfi. (12) 

Thomson scattering complicates the situation by its explicit angle coupling but 
the solution can be obtained with the standard Feautrier scheme. Assuming 
complete frequency-redistribution in spectral lines (Mihalas 1978, p. 29), no ex- 
plicit frequency coupling occurs so that a parallel solution for all frequencies 
enables very efficient vectorization on the computer. 

The following boundary conditions are used for the transfer equation. At 
the inner boundary where the optical depth is at maximum, r = Tmax, we have 

^^Z^) = - Uiyf^irmax), (13) 

Tmax 
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where we specify /J^ = Tmax) from the diffusion approximation: 

1+ -B (14) 

Here B,y is the Planck function and H the nominal (frequency-integrated) Ed- 
dington flux: 

n = <TRT^s/47r, (15) 

with the Stefan-Boltzmann constant ctr. At the outer boundary we take Ti, = 

Tallin = miXjy/2/9, assuming that x is a linear function of m for m < mi. Since 
''"min 7^ 0, it is not exactly valid to assume no incident radiation at the stellar 
surface. Instead we specify = /:/(— /x,rmin) after Scharmer &; Nordlund 
(1982): 

l;;^ = Su{Tmm)[^ - CXp(-rmin//^)] (16) 

which follows from Eq. 8 assuming S{t) = >S'(Tmin) for r < Tmin- Then we get: 

^ / Tjnin 

The boundary conditions are discretized performing Taylor expansions which 
yield second-order accuracy (Mihalas 1978, p. 155). 

3.2. Statistical equilibrium 

The statistical equilibrium equations are set up according to Mihalas (1978, 
p. 127). The number of atomic levels, ionization stages, and chemical species, as 
well as all radiative and collisional transitions are taken from the input model 
atom supplied by the user. Ionization into excited states of the next ioniza- 
tion stage is taken into account. Dielectronic recombination and autoionization 
processes can also be included in the model atom. 

Rate equations As usual, the atomic energy levels arc ordered sequentially by 
increasing excitation energy, starting with the lowest ionization stage. Then for 
each atomic level i of any ionization stage of any element, the rate equation 
describes the equilibrium of rates into and rates out of this level: 

niY.Pij-12'^jPji = 0- (18) 

The rate coefficients Pij have radiative and collisional components: Pij = Rij + 
Cij. Radiative upward and downward rates are respectively given by: 

"''^"^J.diy, (19) 



(20) 
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Photon cross-sections are denoted by CFij{i^)- {ni/nj)* is the Boltzmann LTE 
population ratio in the case of Hne transitions: gi/gj exp(— /iz/jj//cT), where the 
gij are the statistical weights. The LTE population number of a particular level 
is defined relative to the ground state of the next ion, so that in the case of 
recombination from a ground state nf we have by definition (ni/nj)* = ne(f>i{T) 
with the Saha-Boltzmann factor 

UT) = 2.07 • io-^e9irj.-3/2^hu,/kT^ (21) 
9i 

where hvi is the ionization potential of the level i. Care must be taken in 
the case of recombination from an excited level into the next lower ion. Then 
(ni/nj)* = ne4>i ■ 4>i /4>j- 

Dielectronic recombination is included following Mihalas & Hummer (1973). 
Assuming now that j is a ground state of ion k, then the recombination rate 
into level i of ion k — 1 via an autoionization level c (with ionization potential 
/iz^c, having a negative value when lying above the ionization limit) is: 

The reverse process, the autoionization rate, is given by: 

47r2e2 1 

Rij — — ficJ • (2^) 

hmc Vc 

The oscillator strength for the stabilizing transition (i.e. transition i ^ c) is 
denoted by /jc, and J is the mean intensity averaged over the line profile. Our 
program simply takes from the continuum frequency point closest to the 
transition frequency, which is reasonable because the autoionization line profiles 
are extremely broad. The population of autoionization levels is assumed to be 
in LTE and therefore such levels do not appear explicitly in the rate equations. 

The computation of coUisional rates is generally dependent on the specific 
ion or even transition. Several options, covering the most important cases, may 
be chosen by the user. 

Abundance- definition equation The rate equation for the highest level of a 
given chemical species is redundant. It is replaced by the abundance definition 
equation. Summation over all levels usually includes not only NLTE levels but 
also levels which are treated in LTE, according to the specification in the model 
atom. Denoting the number of ionization stages of species k with NION(k), the 
number of NLTE and LTE levels per ion with NL(l) and LTE(l), respectively, 
we can write: 



NION{k) 

E 

1=1 



'NL(l) LTE{1) 
i=l i=l 



n 



kit 



yk(N-ne). (24) 



Uk is the number fraction of element k. 



Model Photospheres with Accelerated Lambda Iteration 



9 



Charge conservation We close the system of statistical equilibrium equations 
by demanding charge conservation. We denote the total number of chemical 
species with NATOM, the charge of ion I with q{l) (in units of the electron 
charge) and write: 



NATOM NION{k) 

E E ^(0 

k=l 1=1 



NUl) LTE{1) 

E E "^kli 

1=1 1=1 



= rie. (25) 



Complete statistical equilibrium equations Wc introduce a vector comprising 
the occupation numbers of all NLTE levels, n = (ni, . . . jHnl)- Then the sta- 
tistical equilibrium equation is written as: 

An = b. (26) 

The gross structure of the rate matrix A is of block matrix form, because tran- 
sitions between levels occur within one ionization stage or to the ground state 
of the next ion. The structure is complicated by ionizations into excited levels 
and by the abundance definition and charge conservation equations which give 
additional non-zero elements in the corresponding lines of A. 

3.3. Radiative equilibrium 

Radiative equilibrium can be enforced by adjusting the temperature stratifica- 
tion either during the linearization procedure or in between ALI iterations. In 
the former linear combination of two different formulations is used and 

in the latter case the classical Unsold-Lucy temperature correction procedure 
(Lucy 1964) is utilized. The latter is particularly interesting, because it allows 
one to exploit the block form of the rate-coefficient matrix. This fact allows an 
economic block-by-block solution followed by a subsequent Unsold-Lucy tem- 
perature correction step. On the other hand, however, this correction procedure 
may decelerate the global convergence behavior of the ALI iteration. 

The two forms of expressing the radiative equilibrium condition follow from 
the requirement that the energy emitted by a volume element per unit time is 
equal to the absorbed energy per unit time (integral form): 

/ xAS.-J.)diy = 0, (27) 
Jo 

where scattering terms in Xv Sind Si, cancel out. In principle, this formulation is 
equivalent to demanding flux-constancy throughout the atmosphere (differential 
form) 

/■oo Q 

/ ■^{UJu)du-H = 0, (28) 

Jo OTv 

where Ti (Eq. 15) is the nominal flux. Here fn is the variable Eddington factor, 
defined as ^ _^ 

/i/ = ^ I^^Uu^diij Uu^dji, (29) 

which is computed from the Feautrier intensity Uy^ (Eq. 9) after the formal so- 
lution. As discussed e.g. in Hubeny (1988) the differential form is more accurate 
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at large depths, while the integral form behaves numerically superior at small 
depths. Instead of arbitrarily selecting that depth in the atmosphere where we 
switch from one formulation to the other, we use a linear combination of both 
constraint equations, which guarantees a smooth transition with depth, based 
on physical grounds (Hubeny & Lanz 1992, Hamann 1994): 



1 r 
iij Jo 



XAS. - Ju) du + A} ^iUJ.) dv - A}n - Fo = 0. (30) 

Jo C)Ti, 



For details on this equation and on our implementation of the Unsold-Lucy 
procedure see Dreizler (this volume). 

We note that explicit depth coupling is introduced by the differential form 
Eq. 28 through the derivative d/dry even if a purely local A* operator is used. 
Therefore the linearization procedure can no longer be performed independently 
at each depth point and the question at which boundary to start becomes rel- 
evant. Numerical experience shows that it is essential to start at the outer 
boundary and proceed inwards. If a tri-diagonal operator is used, nearest neigh- 
bor depth coupling is introduced anyhow. 

3.4. Hydrostatic equilibrium 

We write the equation for hydrostatic equilibrium as Mihalas (1978, p. 170): 

P is the total pressure comprising gas, radiation and turbulent pressures, so 
that: 

with Boltzmann's constant k and the turbulent velocity i^tmb- The hydrostatic 
equation may either be solved simultaneously with all other equations or sepa- 
rately between iterations. The overall convergence behavior is usually the same 
in both cases. If taken into the linearization scheme, and a local A* operator is 
used, then, as in the case of the radiative equilibrium equation, explicit depth- 
coupling enters via the depth derivative d/dm. Again, solution of the linearized 
equations must proceed inwards starting at the outer boundary. The starting 
value at the first depth point (subscript d = 1) is: 

NikTi + ^pit^turbl^T-i) = mi(g-— [ ^^KJykdu] , (33) 
2 \ c Jo pi J 

where hn is the variable Eddington factor denoting the ratio of Hy/Jy at the 
surface, kept fixed during linearization. 

We are also able to account for element separation resulting from pressure 
diffusion, which is the process that governs the spectroscopic appearance e.g. of 
white dwarfs. The numerical method and some results are described by Dreizler 
k Wolff (1999) and Schuh etal. (this volume). 
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3.5. Particle-conservation and fictitious massive-particle density 

The total particle density N is the sum of electron density plus the populations 
of all atomic/ionic states, in all LTE and NLTE levels: 



NATOM NION{k) 

fc=l 1=1 



'NL(l) 



LTE{1) 



E E ^kii 



i=l 



(34) 



The fictitious massive-particle density nn is introduced for notational conve- 
nience. It is defined by: 



NATOM NION{k) 

riH = E ^'<: E 
k=l 1=1 



'NL(l) 



LTE(l) 



E ^kii+ E ^kii 



i=l 



1=1 



(35) 



The mass of a chemical species in AMU is denoted by m^. Introducing the mass 
of a hydrogen atom uih, we may write the material density simply as 



p = nnmH- 



(36) 



3.6. Opacity and emissivity 

Thermal opacity and emissivity arc made up by radiative bound-bound, bound- 
free and free-free transitions. For each species we compute and sum: 



NION 

= E 



'NL{1) NL(l) / \ 



(37) 



1=1 j>i \ 
NL{1) NL{1+1) 

+ E E 'yii^i+iM{nu-nie-'^''/'^^) 

1=1 j>i 

(NL{1+1) LTE(l+l) 
E rn+i^i+ E "HI, 
i=l i=l 

where the total opacity includes Thomson scattering, i.e. Xi/ = Ki/ + neCTg, and 



NION [NL(l) NL{1) 

E E E c.a^..(-)n/.f^e-M"^--.)/'=T 

1=1 I 1=1 j>i 

NL{1) NL{1+1) 

+ E E ^H^/+i,fc(^)nr,e-Wfc^ 

1=1 j>i 

(NL{1+1) LTE{1+1) 
E E 
i=i i=i 



(38) 



(Ta^i+i,fe(z^) denotes the cross-section for photoionization from level i of ion I 
into level k of ion I + I. The double summation over the bound- free continua 
takes into account the possibility that a particular level may be ionized into 
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more than one level of the next high ion. Again, note the definition of the LTE 
population number in this case, which depends on the level {I + l,k) of the 
parent ion: 

* , 4'l+l,l /on\ 

<Pl+\,k 

Note also that the concept of LTE levels (whose occupation numbers enter, 
e.g. the number- and charge-conservation equations) in the atomic models of 
complex ions is therefore not unambiguous. Our code always assumes that LTE 
levels in the model atoms are populated in LTE with respect to the ground state 
of the upper ion. 

The source function used for the approximate radiation transfer is rjy/Ki,, 
thus it excludes Thomson scattering. For the exact formal solution of course, 
the total opacity Xv in the expression Eq. 11 includes the Thomson term (ngiTe). 

3.7. Atomic level dissolution by plasma perturbations 

As high-lying atomic levels are strongly perturbed by other charged particles in 
the plasma they are broadened and finally dissolved. This effect is observable by 
line merging at series limits and must be accounted for in line profile analyses. 
Moreover, line-overlap couples the radiation field in many lines and flux-blocking 
can strongly affect the global atmospheric structure. Numerically, we treat the 
level dissolution in terms of occupation probabilities, which for LTE plasmas 
can be defined as the ratio of the level populations to those in absence of pertur- 
bations. A phenomenological theory for these quantities was given in Hummer 
& Mihalas (1988). The non-trivial generalization to NLTE plasmas was made 
by Hubeny etal. (1994). In practice an individual occupation-probability factor 
(depending on T, ng, and principal quantum number), is applied to each atomic 
level, which describes the probability that the level is dissolved. Furthermore, 
the rate equations Eq. 18 must be generalized in a unique and unambiguous 
manner. For details see Hubeny etal. (1994). 



4. The Accelerated Lambda Iteration (ALI) 

In all constraint equations described above the mean intensities Jy are replaced 
by the approximate radiation field Eq. 2 in order to eliminate these variables 
from the solution vector Eq. 1. In principle the ALO may be of arbitrary form 
as long as the iteration procedure converges. In practice, however, an optimum 
choice is desired in order to achieve convergence with a minimum amount of 
iteration steps. The history of the ALOs is interesting, and was summarized in 
detail by Hubeny (1992). Of utmost importance were two papers by Olson and 
collaborators (Olson et al. 1986, Olson &: Kunasz 1987) who overcame the major 
drawback of early ALOs, namely the occurrence of free parameters controlling 
the convergence process, and thus found the optimum choice of ALOs. Our 
model atmosphere program permits the use of either a diagonal or a tri-diagonal 
ALO, both are set up following Olson & Kunasz (1987). 
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4.1. Diagonal (local) lambda operators 

In this case the mean intensity Jd at a particular depth d in the current iteration 
step is computed solely from the local source function Sd and a correction term 
AJrf, the latter depending on the source functions (at all depths) from the pre- 
vious iteration. Dropping the iteration count and introducing indices denoting 
depth points we can rewrite Eq. 2 as: 

Jd = AldSd + ^Jd- (40) 

In discrete form we can think of A* as a matrix acting on a vector whose elements 
comprise the source functions of all depths. Then ^ is the diagonal element of 
the A* matrix corresponding to depth point d. Writing A^^ = Bd (for numerical 
computation see Eq. 44 below) wc have a purely local expression for the mean 
intensity: 

Jd = BdSd + A J<i. (41) 

4.2. Tridiagonal (non-local) lambda operators 

Much better convergence is obtained if the mean intensity is computed not only 
from the local source function, but also from the source function of the neigh- 
boring depths points. Then the matrix representation of A* is of tri-diagonal 
form and we may write: 

Jd = Cd-iSd-i + BdSd + Ad+iSd+i + AJd , (42) 

where Cd-i and Ad+i represent the upper and lower subdiagonal elements of A*, 
and Sd±i the source functions at the adjacent depths. By analogy the correction 
term becomes: 

AJd = Ad^'Sd' — (Cd-iSd-i + BdSd + Ad+iSd+i). (43) 

Again all quantities for the computation of AJd are known from the previous 
iteration, and the first term gives the exact formal solution of the transfer equa- 
tion. We emphasize again that the source functions in Eq. 42 are to be computed 
from the correct occupation numbers and temperature, all of which are still un- 
known. We therefore have a non-linear set of equations for these quantities, 
which must be solved by cither Ncwton-Raphson iteration or other techniques, 
resulting in the solution of a tri-diagonal linear equation of the form Eq. 5. 

As was shown in Olson et al. (1986) the elements of the optimum A* matrix 
are given by the corresponding elements of the exact A matrix. The diagonal 
and subdiagonal elements are computed following Olson &; Kunasz (1987): 

Ad+i 

1-Bd 

Cd-i 



-ATd_l _ I 



T' 



At, T' 



+ 



1 



-Ard_i 



ATd- 



1] ^ 



(44) 
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with At(/_i = (r^ — Td-ij/n- At large optical depths with increasing Ar steps 
(the depth grid is roughly equidistant in logr) the subdiagonals ^d+i and C^-i 
vanish, and the diagonal approaches unity, reflecting the fact that the radi- 
ation field is more and more determined by local properties of the matter. At 
very small optical depths all elements of A* vanish, reflecting the non-localness 
of the radiation field at these depths. 

4.3. Acceleration of convergence 

Our code takes advantage of an acceleration scheme to speed up convergence 
of the iteration cycle Eq. 2. We implemented the scheme originally proposed 
by Ng (1974) and later by Auer (1987). It extrapolates the correction vector 
(Ji/j^ from the previous three iterations. From our experience the extrapolation 
often yields over-corrections, resulting in alternating convergence or even diver- 
gence. In contrast, use of a tri-diagonal ALO usually results in satisfactorily 
fast convergence, so that the acceleration scheme is rarely used. 

5. Solution of the non-linear equations by iteration 

At each depth the complete set of non-linear equations Eq. 2 for a single it- 
eration step comprises the equations for statistical, radiative, and hydrostatic 
equilibrium, and the particle-conservation equation. For the numerical solution 
we introduce discrete depth and frequency grids. The equations are then lin- 
earized and solved by a suitable iterative scheme. Explicit angle-dependence of 
the radiation field is not required here and is consequently eliminated by the 
use of variable Eddington factors. Angle dependence is considered only in the 
formal solution of the transfer equation. Our program requires an input model- 
atmosphere structure as a starting approximation, together with an atomic data 
file, and a frequency grid. Depth and frequency grids are therefore set up in ad- 
vance by separate programs. 

5.1. Discretization 

After a depth-grid is chosen, we start by computing a gray LTE continuum 
model using the Unsold-Lucy temperature-correction procedure. Depth points 
(typical number is 90) are chosen in equidistant steps on a logarithmic (Rosse- 
land) optical depth scale. The converged LTE model (temperature and density 
structure, given on a column mass depth scale) is written to a file that is read by 
PR02. The NLTE code uses the column mass as an independent depth variable. 

The frequency grid is determined by the atomic-data input file. Frequency 
points are set blue- and red-ward of each absorption edge, and for each spectral 
line. Gaps are filled up by setting continuum points. Finally, the quadrature 
weights are computed. Frequency integrals appearing e.g. in Eq. 30 are replaced 
by quadrature sums, and depth derivatives by difference quotients. 

5.2. Linearization 

All variables x are replaced by x x + 5x where 5x denotes a small perturbation 
of X. Terms not linear in these perturbations are neglected. The perturbations 



Model Photospheres with Accelerated Lambda Iteration 



15 



are expressed by perturbations of the basic variables: 

o o o o NL i~i 

- ox ox . ox ox ^ ^ ox ^ 

As an illustrative example we linearize the equation for radiative equilibrium. 
Most other linearized equations may be found in Werner (1986). Assigning two 
indices {d for depth and i for frequency of a grid with NF points) to the variables 
and denoting the quadrature weights with Wi Eq. 30 becomes: 

NF NF ^ 

'^Wi{^[6Sdi - dJdi] + dxdi[Sdi - Jdi]) + ^^"^ -^{SJdifdi - ^Jd-i,ifd-i,i) 
i=i '^J i=i ^■^^ 

NF NF 

= Fo + A^jH -Y.^i^{Sdi - Jdi) - KY.^^fdiJd^ - fd-uJd-i,i). (46) 

Note that wc do not linearize Atj. Because of this, convergence properties may 
deteriorate significantly in some cases. Perturbations 6Sdi, Sxdi sxe expressed by 
Eq. 45, and the perturbation of the mean intensity Jdi is, according to Eq. 42, 
given through the perturbations of the source function at the current and the 
two adjacent depths: 

5 Jdi = Cd-i,iSSd-i,i + BdiSSdi + Ad+i,iSSd+i,i , (47) 

where A, B, C are the A matrix elements from Eq. 44. The 6Jd-i,i contain 
the term Cd-2,iSSd-2,i which is neglected because we want to account only for 
nearest neighbor coupling. We write 6Sdi with the help of Eq. 45 and observe 
that for any variable z: 

dSdi ^ J_ r drjdj^ _ ^ dx^^ 
dz Xdi \ dz dz J ' 

Derivatives of opacity and emissivity with respect to temperature, electron and 
population densities are computed from analytical expressions (see e.g. Mihalas 
etal. 1975, Werner 1987). We finally get from Eq.46: 



'NF 



Wi dSd-i,i 

d-l,i I 2^ Qrp Xdi^d-l,i 

NF ^ ^ 

+-^jy^ -r^{fdiCd-i,i — fd-i,iBd-i,i) — p^rp \ + 



(NF 



Wi 



dSdi n X , dx, 



Qrp Xdi(l - Bdi) + -^{Sdi - Jdi) 

+ 



+^J^-^{fdiBdi - fd-l,iAdi)-Q^^ 
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\ i 

'^^J -^,{fdiAd+i,i - fd-i,iBd+i,i)—^r^ \ + 

NL NL NL 

1=1 1=1 1=1 

= r.h.s. (49) 

The curly brackets {• • •} denote terms that are similar to those multiplied by 
perturbations of the temperature. Instead of partial derivatives with respect to 
T, they contain derivatives with respect to rie and the populations n^. They all 
represent coefficients of the matrices a, (3, 7 in Eq. 5. 

5.3. Newton-Raphson iteration 

As described in Sect. 2. the linearized equations have a tri-diagonal block-matrix 
form, see Eq. 5. Inversion of the grand matrix (= T sized {NN ■ ND) x {NN ■ 
ND), i.e. about 10^ x 10^ in typical applications) is performed with a block- 
Gaussian elimination scheme, which means that our iteration of the non-linear 
equations is a multi-dimensional Newton-Raphson method. The problem is 
structurally simplified when explicit-depth coupling is avoided by the use of 
a local ALO, however, the numerical effort is not reduced much, because in 
both cases the main effort lies with the inversion of matrices sized NN x NN. 
The Newton-Raphson iteration involves two numerically expensive steps, first 
setting up the Jacobian (comprising a,/9,7) and then inverting it. Addition- 
ally, the matrix inversions in Eq. 6 limit their size to about NN = 250 because 
otherwise numerical accuracy is lost. Two variants recently introduced in stellar 
atmosphere calculations are able to improve both numerical accuracy and, best 
of all, computational speed. 

5.4. Alternative fast solution techniques for non-linear equations: 
Broyden— and Kantorovich— variants 

Broyden's (1965) variant belongs to the family of quasi- Newton methods; it 
was first used in model-atmosphere calculations in Dreizler & Werner (1991), 
Hamann et al. 1991, Koesterke et al. (1992). It avoids the repeated set-up of the 
Jacobian by the use of an update formula. In addition, it also gives an update 
formula for the inverse Jacobian. In the case of a local ALO the solution of the 
linearized system at any depth is: 

S^k = /3k^ Ck. (50) 

Let be the k-th iterate of the inverse Jacobian, then an update can be found 
from: 

1 _ 1 (sfc - /3fc Vfc)0(s^/3fc ^) 

f-'k+l — Pk ^ r,o-l ' 
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where (8> denotes the dyadic product, and where we have defined: 

Sfc = Stpj^ solution vector of preceding hnearization, 

Yfe = Cfe+i — Cfe difference of actual and preceding residuum. 

The convergence rate is super-linear, i.e. slower than the quadratic rate of the 
Newton-Raphson method; but this is more than compensated by the tremen- 
dous speed-up for a single iteration step. It is not always necessary to begin the 
iteration with the calculation of an exact Jacobian and its inversion. Experience 
shows that in an advanced stage of the overall (ALI-) iteration Eq. 2 (i.e. when 
corrections become small, of the order 1%) we can start the linearization cycle 
Eq. 51 by using the inverse Jacobian from the previous overall iteration. Com- 
putational speed-up is extreme in this case, however, it requires storage of the 
Jacobians of all depths. 

More difficult is its application to the tri-diagonal ALO case. Here we 
have to update the grand matrix T which, as already mentioned, is of block 
tri-diagonal form. We cannot update its inverse, because it is never computed 
explicitly. Furthermore we need an update formula that preserves the block tri- 
diagonal form which is a prerequisite for its inversion by the Feautrier scheme 
Eq. 6. Such a formula was found by Schubert (1970): 

T,,, = T, + i^:^^^^^, (52) 

where = Zs^ with the structure matrix Z as defined by: 

_ r 1 if Tij^o 

~ \ if Tij = 0. 

The vectors and are defined as above but now they span quantities over 
all depth instead of a single depth point. With this formula we obtain new 
submatrices a, (3, 7 and c with which the Feautrier scheme Eq. 6 is solved again. 
This procedure saves the computation of derivatives. Another feature realized 
in our program also saves the repeated inversion of q = (/3^ -|- •y^'Dd^i) by 
updating its inverse with the Broyden formula Eq. 51. Similar to the diagonal 
ALO case it is also possible to pass starting matrices from one overall iteration 
Eq. 2 to the next for the update of T and the matrix q""*". In both cases the user 
specifies two threshold values for the maximum relative correction in 5ip which 
cause the program to switch from Newton-Raphson to Broyden stages 1 and 2. 
During stage 1, each new overall cycle Eq. 2 is started with an exact calculation 
and inversion of all matrices involved, and in stage 2 these matrices are passed 
through each iteration. 

Another variant, the Kantorovich (1949) method has been introduced into 
model-atmosphere calculations (Hubeny & Lanz 1992). It is more simple and 
straightforward to implement. This method simply keeps the Jacobian fixed 
during the linearization cycle; it is surprisingly stable. In fact it turns out to be 
even more stable (i.e. it can be utilized in an earlier stage of iteration) than the 
Broyden method in the tri-diagonal ALO case. The user of PR02 may choose 
this variant in two stages in analogy to the Broyden variant. It was found that 
in stage 2 it is necessary to update the Jacobian every 5 to 10 overall iterations 
in order to prevent divergence. 
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6. NLTE metal-line blanketing 

Despite the increase in capacity for solving NLTE model-atmosphere problems 
given by the ALT method combined with pre-conditioning techniques, the blan- 
keting by millions of lines from the iron-group elements arising from transitions 
among some 10^ levels can be attacked only with statistical methods. These 
were introduced into NLTE model-atmosphere work by Anderson (1989, 1991). 
At the outset, model atoms are constructed by combining many thousands of 
levels into a relatively small number of superlevels which can be treated by ALI 
(or other) methods. Then, in order to reduce the computational effort, two ap- 
proaches were developed which vastly decrease the number of frequency points 
(hence the number of transfer equations to be solved) needed to describe prop- 
erly the complex frequency dependence of the opacity. These two approaches 
have their roots in LTE modeling techniques, where statistical methods are ap- 
plied for the same reason in the treatment of opacity: The Opacity Distribution 
Function (ODF) and Opacity Sampling (OS) approaches. Both are based on 
the fact that the opacity (in the LTE approximation) is a function of two only 
local thermodynamic quantities. Roughly speaking, each opacity source can be 
written in terms of a population density and a photon cross-section for the re- 
spective radiative transition: ~ niaiuiy)- In LTE the population follows from 
the Saha-Boltzmann equations, hence n; = ni{ne,T). The OS and ODF meth- 
ods use such pre-tabulated (on a very fine frequency mesh) Kj/(ne, T) during the 
model atmosphere calculations. The NLTE situation is more complicated, be- 
cause pre-tabulation of opacities is not useful. The occupation numbers at any 
depth now also depend explicitly on the radiation field (via the rate equations 
which replace the TE Saha-Boltzmann statistics) and thus on the populations 
in every other depth of the atmosphere. As a consequence, the OS and ODF 
methods are not applied to opacity tabulations, but to tabulations of the photon 
cross-sections cj(i^). These do depend on local quantities only, e.g. line broad- 
ening by Stark- and Doppler-effects is calculated from T and rig. In the NLTE 
case cross-sections take over the role which the opacity played in the LTE case. 
So, strictly speaking, the designation OS and ODF is not quite correct in the 
NLTE context. 

The strategy in our code is the following. Before any model atmosphere 
calculation is started, the atomic data are prepared by constructing superlevels, 
and the cross-sections for superlines. Then these cross-sections are either sam- 
pled on a coarse frequency grid or ODFs are constructed. These data are put 
into the model atom which is read by the code. The code does not know if OS 
or ODFs arc used, i.e. it is written to be independent of these approaches. 

Further details about our implementation are described in Dreizler & Wer- 
ner (1993), Haas etal. (1996), Werner & Dreizler (1999), and by Rauch & Deet- 
jen in this volume. 



7. Summary 

We have presented in detail our technique for numerical solution of the classical 

model-atmosphere problem. The construction of metal-line-blanketed models in 
hydrostatic and radiative equilibrium under NLTE conditions was the last and 
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longest-standing problem of classical model atmosphere theory. It has finally 
been solved with a high degree of sophistication. The essential milestones for 
this development, starting from the pioneering work of Auer &: Mihalas (1969) 
are: 

• Introduction of Accelerated Lambda Iteration (ALI, or "operator split- 
ting") methods, based upon early work by Cannon (1973) and Scharmer 
(1981). The first ALI model atmospheres were constructed by Werner 
(1986). 

• Introduction of statistical approaches to treat the iron-group elements in 
NLTE by Anderson (1989, 1991). 

• Linear formulation of the statistical-equilibrium equations (Rybicki &: Hum- 
mer 1991, Hauschildt 1993). 

• Computation of atomic data by Kurucz (1991), by the Opacity Project 
(Seaton et al. 1994) and subsequent improvements, and by the Iron Project 
(Hummer etal. 1993). 
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